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ON THE CONSTRUCTION OF VIRTUAL INTERIOR POINT 
SOURCE TRAVEL TIME DISTANCES FROM THE HYPERBOLIC 
NEUMANN-TO-DIRICHLET MAP 


MAARTEN. DE HOOP PAUL KEPLEY *, AND LAURI OKSANEN § 

Abstract. We introduce a new algorithm to construct travel time distances between a point 
in the interior of a Riemannian manifold and points on the boundary of the manifold, and describe 
a numerical implementation of the algorithm. It is known that the travel time distances for all 
interior points determine the Riemannian manifold in a stable manner. We do not assume that there 
are sources or receivers in the interior, and use the hyperbolic Neumann-to-Dirichlet map, or its 
restriction, as our data. Our algorithm is a variant of the Boundary Control method, and to our 
knowledge, this is the first numerical implementation of the method in a geometric setting. 
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1. Introduction. We consider the inverse boundary value problem for the acous¬ 
tic wave equation on a domain or on a Riemannian manifold M. The problem is to 
recover a spatially varying sound speed c inside the domain given the corresponding 
Neumann-to-Dirichlet map, or its restriction on a part of the boundary of the domain, 
say r c dM. 

In the acoustic case, waves propagate in a domain M with speed c, and the travel 
time of a wave between a pair of points x and ?/ in M is given by the Riemannian 
distance d{x,y) when computed with respect to the travel time metric c~^dx^. For 
this reason, it is natural to formulate the inverse boundary value problem by using 
concepts from Riemannian geometry. This also allows us to consider anisotropic sound 
speeds in a unified way. 

We present a new method to use the local restriction of the Neumann-to-Dirichlet 
map to determine travel time distances of the form d{x, y) where a; is a point in the 
interior of the domain and y S F, that is, y is a point in the set where we have 
boundary measurements. We refer to these travel times as point source travel time 
data, since the distance d(x, y) corresponds to the first arrival travel time from a 
(virtual) interior point source located at x as recorded at the boundary at y. We 
emphasize that our method synthesizes the travel times from a point source in the 
interior of M without requiring an actual receiver or source at that location. 

Using the travel time distances we can introduce Dirac boundary sources which 
generate waves the singularities of which focus in a point. Such boundary sources 
have been referred to as focusing functions in reflection seismology m- 

To motivate our results, we note that travel times have previously been shown 
to determine a Riemannian manifold with boundary see e.g. m- In particular, in 
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the full boundary data case, it has been shown that this determination is even stable 
m- Furthermore, it can be shown that travel times determine shape operators that 
appear as data for the generalized Dix method mi. This is of particular interest in 
the isotropic case, since that method allows for the local nonlinear reconstruction of 
a sound-speed near geodesic rays. 

Our method to determine boundary distances works by first using a variant of 
the Boundary Control method (BC method) to determine volumes of subdomains of 
M referred to as wave caps. The volumes are computed by solving a collection of 
regularized ill-posed linear problems. The BC method originates from [^, where it 
was used to solve the inverse boundary problem for the acoustic wave equation. We 
note that [8] was the first variant of the BC method posed in a geometric setting. 
For a thorough overview of the BC method, we refer to [H] and [7]. Regularization 
theory was first combined with the BC method in m, and the procedure that we use 
to determine volumes was developed in |28) . 

In the present paper, we introduce a procedure to construct point source travel 
time data from the volumes of wave caps, and hence reduce the inverse boundary 
value problem to the stable problem of determining the Riemannian manifold (M, g) 
given the point source travel time data. In particular, our procedure splits the inverse 
boundary value problem to an ill-posed but linear step (the volume computation) 
and a non-linear but well-posed step (distance estimation and reconstruction of the 
manifolds). 

We describe a computational implementation of our method, and provide a nu¬ 
merical example to demonstrate our technique. We remark that our numerical ex¬ 
ample provides the first computational realization of a geometric variant of the BC 
method. Moreover, we explain how the instability of the volume computation step is 
manifest in our numerical examples. All variants of the BC method that use partial 
data contain similar unstable steps, and we believe that the instability of the method 
reflects the ill-posedness of the inverse boundary value problem itself. However, con¬ 
trary to Calderon’s problem, that is, the elliptic inverse boundary value problem 
[HIM], it is an open question whether the inverse boundary value problem for the 
wave equation is ill-posed in general. For results concerning the stability of Calderon’s 
problem we refer to [U US] . 

Also, we point out that under favorable geometric assumptions, the hyperbolic 
inverse boundary value problem is known to have better stability properties than 
Calderon’s problem, see e.g. [SllMlEllEaisl] and references therein. However, these 
geometric assumptions do not apply to the partial data problem that we consider 
here. 

2. Statement of the Results. In this section we describe our data and as¬ 
sumptions, and what we intend to recover from the data. 

2.1. Direct problem as a model for measurements. We work in the fol¬ 
lowing setting. 

Assumption 1. M is a smooth compact connected manifold with smooth bound¬ 
ary DM, and g is an unknown smooth Riemannian metric on M. 

Let fi G C°°(M) be a strictly positive weight function, note that we do not 
require any additional information about g. We consider the following wave equation 
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on {M,g) with boundary source / S C“((0,oo) x dM), 

dfu{t,x) — Ag^f^u{t,x) = 0, {t,x) G (0,oo) X M, 

(2.1) Ngu{t,x) = f{t,x), (t, x) G (0, oo) X r 

m( 0, •) = 0, i9tu(0, •)=0, X G M. 

The operator is the weighted Laplace-Beltrami operator, given by, 

Ag^f,w{x) := divg(/xgradg(u;)). 


and, Ng, is the associated Neumann derivative, 


NgW := -g, {v, gradg(w))g , 


where ly is the inward pointing unit normal vector to dM in the metric g, divg and 
gradg are respectively the divergence and gradient on {M,g), and (•, ■)g denotes the 
inner product with respect to the metric g. We remark that Ag^g is defined so that 
we can also handle the usual acoustic wave equation. Indeed, if we consider a domain 
M C M" along with a strictly positive function c G C°°{M), then with respect to the 
conformally Euclidean metric g = c~'^dx'^, the weight g = yields Ag g^ = c^A. 


We now introduce our data model. We denote the solution to (2.11 with Neumann 


boundary source / by u^{t,x). For T > 0, and T C dM an open set, we define the 


local Neumann-to-Dirichlet operator associated with (2.1) by, 


^ ^^■^l(o,2T)xrj / G Co°((0:2T) X r). 


The map extends to a bounded operator on L^((0, 2T) x T), see for instance 
|22| . For data, we suppose that Ap^ is known. We note that Ap^ models boundary 
measurements for waves generated with acoustic sources and receivers located on 
T, where the waves are both generated and recorded on T for 2T units of time. In 
addition, we note that in seismic applications the Neumman-to-Dirichlet map appears 
in simultaneous source acquisition. 

Our primary interest is in constructing distances with respect to the Riemannian 
metric g, and we denote the Riemannian distance between the points x,y G M hy 
d{x,y). To simplify our distance computation procedure, we assume: 

Assumption 2. The distances d{y, z) are known for y,z GT with d{y, z) < T. 
We note that this is not a major limitation since for y, z G T with d{y, z) < T, the 
data Ap"^ determines d{y,z), see e.g. [T51 Section 2.2]. 

2.2. Reconstruction of the point source travel time data. We define 
R{M) to be the set of boundary distance functions on M, 


(2.2) R{M) = {rj; : x G M, and for z G dM, r^iz) := d{x, z)}. 


We note that, for x G M and z G dM, r^iz) gives the minimum travel time from x 
to z. With this interpretation, rx(z) represents the first arrival time at z from a wave 
generated by a point source located at x. 

In Section]^ we develop a method to synthesize values of from Ap^ for points x 
indexed by a set of coordinates known as semi-geodesic coordinates^ We refer to this 
procedure as forming point source travel time data, since our procedure reproduces 

^Such coordinates are considered in seismology, where they are referred to as image ray coordi¬ 
nates m 
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the travel time information for a point source located at x without having a source 
or receiver there. 

The geometry of the supports of solutions to (2.1) inform our constructions. To 
be explicit, let r be a function on T satisfying 0 < r(z) < T for all 2 € T and define 
the domain of influence of r, 


M(t) := {x S M : there exists a z G T such that d{x,z) < t(z)} . 


We depict M(t) in Figure la 


Consider the set 


(2.3) := {{f z) e (0, T) X f : t & {T - t{z), T)}. 


We recall that solutions to in exhibit finite speed of propagation in the metric g, 
and specifically, if supp(/) C Sr then supp(M'^(r, •)) C M{t). 

When r is a multiple of an indicator function, we will occasionally use a special 
notation for M(t). To be specific, we denote the indicator function of a set S' by Is, 
and for s > 0 we will use the notation M(r,s) := M(slr) and M{y,s) := M(sl{y}). 

We denote the unit sphere bundle SM := G TM : |^|g = 1}, and define the 
inward/outward pointing sphere bundles by d±SM := G dSM : (^,±^)g > 0}, 
where v is the inner unit normal vector field on dM. We define the exit time for 
(x, f) G SM \ d+SM, by 


tm(x, i) := inf{s G (0, oo) : 7(5; x, G dM), 


where 7 (-;x,^) is the geodesic with the initial data 7 ( 0 ) = x, 7 ( 0 ) = 

For y G F we define crp(y) to be the maximal arc length for which the normal 
geodesic beginning at y minimizes the distance to F. That is. 


(Tr(y) := max{s G ( 0 ,TM(y, i^)] : d{"f{s-,y,iy),T) = s}. 

We recall, see e.g. [HI p. 50] that crr(y) > 0 for y G F. Moreover, or is lower 
semi-continuous, see e.g. [531 Lemma 12]. We define 

(2.4) x(y, s) := 7(5; y, iz) for y G F and 0 < s < crr(y)- 

The mapping (y, s) 1— >■ x(y, s) is a diffeomorphism from {(y, s) : y G F, 0 < s < crr(y)} 
onto its image, so we will refer to (y, s) as semi-geodesic coordinates for x(y, s). This 
is a slight abuse of terminology, since the pair (y, s) belongs to F x [ 0 ,oo) instead 
of a subset of K". On the other hand, by selecting local coordinates on F these 
“coordinates” can be made into legitimate coordinates. 

Next, we recall the definition of the cut locus of F, 


C = {x(y,crr(y)) : y G F} . 

We depict C in Figure]^ Due to the lower semi-continuity of tTp and the boundedness 
of F, one sees that the distance between F and C is positive. 

We will use the following notions of volume: let dVg and dSg denote the Rieman- 
nian volume densities of {M,g) and {dM, g\gM) respectively. We remark that dSg is 
determined on F by Ap^, see e.g. m Section 2.2] so we assume that it is known. We 
define the natural Riemannian volume density associated with Ag g by dVg := gdVg. 
We remark that its name derives from the fact that ^g,g is self-adjoint on Lf{M-, dVg) 
with domain iLg (M)niL^(M). The volume density dVg determines a volume measure 
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Fig. 1: (a) The domain of influence for a r in C(r) along with the profile of r. (b) The 
cut locus of r along with a pair of equal length geodesics showing the break-down 
of the semi-geodesic coordinates at C. The shaded region is the subset of M that 
supports semi-geodesic coordinates. 


which we denote Vol^. In addition, we will use the following shorthand notation for 
volumes of domains of influence mir) := Vol^(M(r)). 

We now describe a set of geometrically relevant subsets whose volumes will allow 
us to determine distances. Let y G T, and let s,h > 0 satisfy s + h < crr(y)- We define 
the wave cap, 

capr(y, S,h) := M{y,s + h) \ M°(r,s), 

where M°(r, s) = {x G M : d(x,r) < s}. Note that under the above hypothe¬ 
ses, x{y,s) belongs to capp(y, s, h). We will use volumes of wave caps to determine 
distances. 

Our main result is an algorithm to use the data Ap^ to construct distances of 
the form r,^^y^s){z) for y,z G T and s > 0 with d{x(y,s),z) < min(crp(y),T). Our 
procedure can also be viewed as a constructive proof of the following known result, 
see e.g. [T^ : 

Theorem 2.1. Let y, z G T and s > 0 with d{x{y,s),z) < min((Tp(y),T). Then 
AP’ determines rx{y,s)(z)- 

The constructive proof will be given in Section We note that this construction 
can also be viewed as a series of experiments. Following the proofs in Section [T2l we 
provide an algorithmic overview of our distance computation procedure. 

3. The Boundary Control method. In this section we describe the elements 
of the BC method required to determine m(r) from Ap. In addition, we briefly 
contrast our technique to alternative approaches to the BC method, and provide an 
overview of some computational aspects of the BC method. 

The purpose of the BC method is to gain information about the interior of M 
by processing boundary measurements for waves that propagate in M. To begin, 
we recall the control map, Wt, which takes a Neumann boundary source / to the 
corresponding solution at time T. That is, let r : F —>• [0, T] be continuous or a step 
function with open level sets and define, 

Wrf := uf{T, •), Wr : L^(Sr) -G L^(M). 

The map Wr is continuous, compact even, as a map from L‘^{Sr) into L?{M), see e.g. 
[221 13S]- We will write W := Wr in the special case that t = T, and we note that in 
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this case, Sr = (0,T) x F. Thus, for any other r : F —>• [0,T], Wr can be viewed as a 
restriction of W to sources supported in Sr ■ 

One cannot directly observe the output of Wr from boundary measurements be¬ 
cause its output is a wave in the interior of M. Thus, in order to deduce information 
about the interior of M, one forms the connecting operator, 

Kr := WrWr, K : L^{Sr) L^{Sr). 


The continuity of Wr implies that Kr is a continuous operator on Lp'{Sr)- The practi¬ 
cal utility of Kr is that it can be computed by processing the boundary data, Ap^, see 


(3.2) below. This fact was first observed by Blagoveschenskii in the l-|-l-dimensional 
m- We remark that Kr derives its name from the fact that it connects the 


case 


inner product on boundary sources with the inner product on waves in the interior. 
That is, for /, h in C^{Sr), 


(3.1) 


(u-^(T, •),M^(T, •))i2(M;dV,,) = {f,Kh)L2(^Sr-,dt0dSg 


We next recall the “Blagoveschenskii Identity,” which gives an expression for Kr 
in terms of the data Ap^. In particular, we use the expression for Kr appearing in 

m. 

(3.2) Kr = Pr {jA'fe - RA'^Rje) Pr- 

Here, 0 : L^((0,T) x F) —)• L^((0, 2r) x F) is the inclusion (zero padding) given by: 


e/(t,-) 


/(t,-) 0<t<T, 

0 T <t<2T, 


R : L^((0,T) X F) — >■ L‘^{{0,T) x F) is the time reversal on (0,T) given by: 


Rf{t,-):=f{T-t,-) 0<t<T, 

J : L^((0, 2r) X F) —>• L^((0,T) x F) is the time integration, given by: 

1 

f{sr)ds 0<t<T, 

and Pr : T^((0,T) x F) —L^((0,r) x F) is the orthogonal projection onto L'^{Sr) 
given by: 


(3.3) 


Prf ■= ISr ■ /■ 


We will use the special notation K := Kr when t = T. In this case, the operator 
Pr coincides with the identity and (3.1) can be written as AT = JAp^0 — RA^RJQ. 
Thus, for any r, the operator Kr can be expressed as Kr = PrKPr. 


3.1. Overview of BC method variants. There are several variants of the BC 
method, all of which are based on solving control problems of the form: Given a 
function (j) on M, and a function r : F —> [0, T], find a boundary source / such that 


(3.4) 


Wrf = 
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In general, this problem is not solvable since the range of WV is generally not closed. 
On the other hand, it can be shown that approximate controllability holds, that is, 
there is a sequence C CQ°{Sr) such that 


(3.5) 


lim Wrfj = in L^{M). 


The approximate controllability follows from the hyperbolic unique continuation result 
by Tataru ^S] by a duality argument, see e.g. m p-157]. 

The original version of the BC method [5] uses the Gram-Schmidt orthonormal¬ 


ization to find a sequence satisfying (3.51. The method was implemented 


numerically in [5], and it requires choosing an initial system of boundary sources, see 
step 2 in [5l p. 233]. No constructive way to choose the initial boundary sources is 
given, and some choices may lead to an ill-conditioned orthonormalization process, 
see the discussion in m- 

More recently, Bingham, Kurylev, Lassas and Siltanen introduced a variant of the 
BC method where the Gram-Schmidt process is replaced by a quadratic optimization 
m- Their method is posed in the case T = dM, and is based on constructing a 


sequence (/j)^i such that the limit (3.51 becomes focused near a point. To elaborate, 
their method considers an arbitrary h G L^{{0,T) x dM) with </> chosen as (j) = Wh. 
For a point y G dM and small enough 0 < s,r < T, they choose appropriate r 
to produce a sequence of sources C Sr such that Wfj -G ^ca.pgi^{y,s,r)Wh. 

However, no constructive procedure to choose the boundary source h is given, and 
some choices may lead to sequences such that this limit vanishes also near the point 
where it should be focused, see the assumption on the non-vanishing limit in [B 
Gorollary 2]. We note that the method [TU] has not been implemented numerically. 

Our approach employs a quadratic optimization similar to |10j but differs from it 
by selecting (p = 1 in place of Wh. By solving the approximate control problems for 
this choice of (/>, we can compute volumes m(T) for certain functions r : T —> [0, T]. We 
note that the method we use to compute these volumes was developed in [251HH] : and 
it was applied to an inverse obstacle problem in |30j . Here we show how to compute 
the boundary distance functions from the volumes TO(r). 

Our method contains only constructive choices of boundary sources, and it allows 
us to understand the numerical errors that we make in each step of the algorithm. In 
Section we see that the dominating source of error in our numerical examples is 
related to the instability of the control problem (3.4) under the constraint supp(/) C 
Sr ■ This instability is inherent in all the variants of the BG method mentioned above. 

In addition to j^, the only multidimensional implementation of a variant of the 
BG method, that we are aware of, is m- This variant is based on solving the control 
problem (3.4) without the constraint supp(/) C Sr- The target function p is chosen to 
be harmonic, and the method exploits the density of products of harmonic functions 
in L^{M). Such an approach works only in the isotropic case, that is, in the case of 
the wave equation d^ — c(x)^A where the sound speed c{x) > 0 is scalar valued. 

We also mention that the original version of the BC method |6] assumes the wave 
equation to be isotropic, and that in [24] . an approach similar to [31] was shown 
to recover a lowpass version of the sound speed in a Lipschitz stable manner under 
additional geometric assumptions. Furthermore, we refer to [18] for a comparison of 
the BC method and other inversion methods in the 1-|-1-dimensional case. 


3.2. Regularized estimates of volumes of domains of influence. We now 

explain how we pose our approximate control problems, and how we use their solutions 
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to compute volumes of domains of influence. To begin, let r : T —?> [0, T] be either a 
step function with open level sets or r S C(r). We obtain an approximate solution 
to (3.4) with right-hand side (/) = 1, by solving the following minimization problem; 


for a > 0 let 
(3.6) 


fa ■= argmin 




L^iM-dV^) 


\L^{S^;dt(»dSg)- 


As was shown in 


for T as above: this problem is solvable, the solution can be 
obtained by solving a linear problem involving Kr, and u^°‘{T, •) —)■ 1m{t) as a —)■ 0. 
For the convenience of the reader, we outline the proof here, and moreover, we recall 
that the approximate control solutions, fa, can be used to compute m(T). 


To show that (3.6) has a solution we first recall two results about Tikhonov 
regularization. For proofs see e.g. pn Th. 2.11] and [20], respectively. 

Lemma 3.1. Suppose that X and Y are Hilbert spaces. Let y G Y and let 
A : X ^ Y be a bounded linear operator. Then for all a > 0 there is a unique 
minimizer of 

\\Ax - yf + a\\xf 


given by Xa = {A*A + a)~^A*y. 

Lemma 3.2. Suppose that X and Y 
A : X ^ Y be a bounded linear operator 
a —>■ 0, where Xa = (A*A + a)~^A*y, a > 
projection. 

Since Wr is bounded, the first Lemma 
the second lemma to our current setting, 
compute Wfl. Toward that end, we recall 
of propagation. When r is a step function 
that the inclusion 


are Hilbert spaces. Let y G Y and let 
with range R{A). Then Axa —t Qy as 
0, and Q : Y ^ ihe orthogonal 


implies that (3.6) is solvable. To apply 
we must describe the range of Wr and 
that supp(lF,-/) C M(t) by finite speed 
Tataru’s unique continuation [OS] implies 


(3.7) 


{Wrf; f G L\Sr)} C L^(M{t)), 


is dense, see e.g. m Th. 3.10]. The result was extended to the case of r G C'(r) in 
[28] . Thus R{Wr) = L'^{M{t)) for t he fu nctions r under consideration. To compute 
1, we note an equality similar to (3.1) that is satisfied for / G Lf‘{Sr)'. 


W: 


(3.8) 


(w^(T, •), l)i2(M;£;v^) = {f,PTb)L^((o,T)xdM-dt(»dSg)- 


Here, Pr is defined by ( ]3.3| ), and b{t,x) := T — t for {t,x) G (0,T) x F. Thus, 

Wfl = Prb. 


Applying Lemmas 3.1 and 3.2 to the observations above, we see that for each 
a > 0, equation (3.6) has a unique solution fa, given by: 


(3.9) 


fa := {WfWr + ay^Wfl = (Kr + a)-^Prb. 


thus fa is obtained from the data. Moreover, the waves Wrfa satisfy Wrfa —t QtI 
in Lf‘{M) as a tends to zero, where Qr is the projection of Lf‘{M) onto the s ubsp ace 
R{Wr) = L‘^{M{t)). Note that = lM(r)- Using this fact and applying (3.8) to 
fa we conclude. 


(3.10) m{T) = lim (/a,T’T&)L2((0,T)x3M;dt(g.dS„)- 

Thus we can compute m{T) from operations performed on the data Ap^. 
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4. Constructing distances. In this section, we present our proof of Theorem 
|2.1[ We accomplish this through a sequence of lemmas that are designed to illuminate 
the steps required to turn the theorem into an algorithm. In addition, we provide 
an alternative technique to determine distances, which we use in our computational 
implementation. 

4.1. Constsructive proof of Theorem |2.1[ The following lemma provides a 
bound on the distance between a point and a wave-cap. 

Lemma 4.1. Let ?/ S T, s G {0,ar{y)), and h G (0,CTr(2/) — s). Let z G T and 
r > 0. Then d{z, capp( 2 /, s, h)) < s + r if and only if 

(4.1) m(slr -I- rlz -I- hly) — m(slr -f rl^) < TO(slr -f hly) — m(slr). 

We note that tests whether there is an overlap between the sets capp(2/, s, h) and 
capp(?/, s, r), see Figure^ 

Proof As h > 0 and h < tTp(j/) — s, we see that capp(?/, s,/i) contains a 
non-empty open set. In particular, it has strictly positive measure. Moreover, if 
d{z, capp(y, s, h)) < s + r then the intersection of capp( 2 /, s, h) and M{z,s + r) con¬ 
tains a non-empty open set and has strictly positive measure. 

Notice that to(s1p -I- hly) is the measure of M{y,s + h) U M(r, s) and that 
to(s1p -I- hly) — m(slp) is the measure of capp(y, s, h). Indeed, 

Yo\y,{M{y,s + h) U M(r, s)) = Yo\y,{{M{y, s + h) U M(r, s)) \ M(r, s)) 

+ Vol^(M(r,s)) 

= Vol^(capp(y, s,/i)) -f Vol;,(M(r, s)). 

Analogously, m(slp + rlz + hly) — m(slp + rlz) is the measure of 

M{y,s + h)\ (M(r, s) U M{z,s + r)) = capp(y, s, h) \ M{z, s + r). 


If d{z, capp(y, s, h)) < s + r then the intersection of capp(y, s, h) and M{z,s + r) has 
strictly positive measure, whence (4.1) holds. 

On the other hand, if d{z, capp(y, s, h)) > s + r then capp(y, s, h) D M{z,s + r) 
is contained in the topological boundary of M(z,s + r) which is of zero measure [55]. 
Thus 


i(slp -I- rlz -I- hly) — m(slp -|- rlz) = to(s1p -I- hly) — m(slp). 


and (4.1) does not hold. 


□ 


The next lemma demonstrates that when s < crp(y), the wave caps capp(y, s, h) 
tend, in a set-theoretic sense, towards x(y, s). 

Lemma 4.2. Let y G T and s G (0,(Tp(y)). Then, 


(4.2) 


n capp(y, s,/i) = {a;(y,s)}. 


h>0 


Proof Let y, s as above, and let /(y, s) denote the left hand side of (4.2). Let w be 
any point belonging to J(y, s). Then w G capp(y, s, h) for all h > 0, so s < d{T, w) and 
d{y, w) < s + hior all h > 0, thus d(y, w) < s. Since y G T, we conclude s = d{T, w) = 
d{y,w). On the other hand, if w is a point in M satisfying s = d{T,w) = d{y,w), 
then w G capp(y, s, h) for any h > 0, hence w G /(y, s). We conclude. 


(4.3) 


J(y, s) = {w G M : d{y,w) = d{r,w) = s}. 
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(b) s + r < ^( 2 , capp(^, s,/i)) 


Fig. 2: A wave cap Fig. 3: The light gray regions indicate the wave caps used in 
Lemma |4.1| and the dark gray region indicates the overlap 
between the caps. 


Because s < ariy), we have d(x{y, s),y) = d{x{y, s),T) = s, so x{y,s) G I(y,s). It 
remains to show that no other points belong to I{y, s). 

Let w belong to I{y,s), we will show that w = x{y,s). If we knew for certain 
that w belonged to the image of the semi-geodesic coordinates, then this would be 
immediate from the definition of these coordinates. On the other hand, if we did not 
require F to be open, then simple examples show that for points y in the topological 
boundary of F it is possible that I{y, s) has many points. We demonstrate that when 
F is open this cannot happen. 

Since M is a compact connected metric space with distance arising from a length 
function, the Hopf-Rinow theorem for length spaces applies and we conclude that there 
is a minimizing path /3 : [0,1] —>■ M from y to w. By [2], /3 is and we may assume 
that it is unit speed parameterized. Hence I = s. As /3 is minimizing from both y and 
F to w, we see that /3(0) = v. Thus /3 coincides with x{y,t) for t < mhi(s, tm{ y, ■ 

But s < crr(y), hence s < TM{y, v). Thus we see that w = j3{s) = x{y, s). □ 

We use the preceding lemma to show that, when h is small, the distance between a 
point z GT and the wave cap capp(?/, s, h) surrounding x{y, s) yields an approximation 
to d{z,x{y, s)). We depict this approximation in Figure]^ 



Fig. 4: Cartoon demonstrating that d{z, capp(y, s, hj) -G d{z, x{y, s)) as h —>■ 0. 

Lemma 4.3. Fory,z € F, and s < ar{y), d{z, {y, s,h)) —>■ d{z,x(y,s)) as 
h^O. 

Proof. Let {hj} C K+ be a sequence for which hj } 0. Then each capp(?/, s, hj) is 
compact, so there exists Wj G capp(y, s, hj) such that d{z,Wj) = d(z, capp(j/, s, hj)). 
Because M is a compact manifold, the sequence {wj} has a convergent subsequence 






































CONSTRUCTING TRAVEL TIME DISTANCES 


11 


{wj,_} converging to a point w. Since the wave caps cap^iv, s,h) nest, the tail of 
{wj^.} belongs to the closed set capr(?/, s, hj^) for each hence w G capp(y, s, h) for 
each h > 0. By the previous lemma, we conclude w = x{y, s). 

Together, continuity of the distance function and the particular choice of the 
Wj^ imply that, d{z, cap^iVi s,hj^)) = d(z,Wj^) —>■ d{z,x{y, s)). Since the wave caps 
nest, the sequence {d{z,Wj)} is monotone non-decreasing, and since it is bounded 
above it has a limit. In particular, any subsequential limit coincides with the limit. 
Thus we conclude that d{z, cap^{y, s,hj)) -G d{z,x{y,s)) as j —)■ oo, and in turn, 
d{z, capp(y, s, h)) -p d{z, x{y, s)) as h ^ 0. □ 


The volumes appearing in Lemma 4.1 cannot be computed directly with the r’s 
appearing in the regularized volume determination. That is, the lemma requires us 
to compute volumes such as m(slp -frl^ + hly), but the function r = sip -l-rl^ + hly 
is equivalent to sip in L^((0,T) x T). As a result, the set L‘^{Sr) will produce waves 
that fill L^(M(slp)) as opposed to the desired set L‘^{M{t)). The problem is that the 
spikes hly and rl^ have supports with dSg measure zero. The remedy is to replace 
the spikes by functions that produce the same domains of influence but have better 
supports. To accomplish this, for y GT and R G [0, oo), we define on T by: 


(4.4) 


Tjf(z) :=R-d{z,y) 


for z GT. 


Note that Ty is continuous. We recall that under Assumption the distances d{y,z) 
for y,z G T with d{y,z) < T are known (or, alternatively, that they have been 
computed in some other fashion from Ap^). Thus under our assumptions the functions 
are known. 

Lemma 4.4. Let y,z GT, s,r,h > 0. We will use the notation f W g to denote 
the function obtained by taking the pointwise maximum of f and g. Then, we have 
the following equalities, 


(4.5) 

(4.6) 

(4.7) 


M(t"+'‘ V r. 


M(r;) = M{rly), 

*+-Vs) = M(/ilj,+rl,+ slp), 

s+r 


V s) = M{s1t + rlz). 


Proof Let x G M{rly), then d{y,x) < r. 


Since Ty(y) = r, we have that d{y, x) < 


Tf{x), hence x G M{Tf). Now let x G M{Ty). Then there is a point z G T for 
which d{x, z) < Tf(z). Applying the definition of Ty, we find r > d(x, z) + dr{y, z) > 
d{x,z) +d{y,z) > d{x,y). Hence x G M{rly). We conclude that M{Ty) = M{rly). 
We demonstrate equality (4.6) and note that (4.7) is proved in an analogous 
_ g-s+h Y ^s+r Y Then x G M{t) just in case d{x,p) < t{p) for 


fashion. Let 
some p G T, 


which happens if and only if d{x,p) is less than Ty~^’^{p), (p) , or s. 


The preceding paragraph implies that this happens just in case x belongs to M((s -|- 
/i)ly), M((s-|-r)li), or M(slp), which happens if and only if cc G M{sly + rlz + slr)- 
□ 

We are finally in a position to prove Theorem |2.1[ 


Proof [Of Theorem 2.1 First, let r and h be positive numbers satisfying s -|- r < 


T and s + h < T. Define functions n = s, T2 = V 


T3 = T, 


V s, and 


_ ~s+h 


Ti = r, 


V T?+’' V s. 


Using the regularized volume determination from equation 
(3.10), we compute the volumes m(ri) for i = 1,...,4. Then, Lemma |4.4| implies 
that m(ri) = to(sIp), m(r2) = m(slp -f hly), mir^) = to(s1p -I- rl^), and m(r4) = 
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TO(slr + hly + rlz), thus we have determined the volumes appearing in (4.1). By 
Lemma 4.1 we can compute d(z, cap-p(y, s, h)) by 


(4.8) 


d(z, capp(j/, s, h)) = s + inf{r : 0 < r < T — s, and (4.1) holds}. 


Finally, by Lemma 4.3 we can compute d{z,x{y, s)) by 


(4.9) 


diz, x{y, s)) = lim d{z, cap^iv, s, h)). 
h—^0 


4.2. Alternative distance estimation method. The method to estimate dis¬ 
tances derived from Theorem |2.1| uses the fact that, under the hypotheses of the the¬ 
orem, the distance between a point z G T and the wave cap capp(j/, s, h) serves as 
an approximation to d{z,x{y,s)), and that this approximation improves as h —>■ 0. 
However, in the case where g is the Euclidean metric, d{z, capp( 2 /, s, h)) converges to 
d{z,x{y,s)) with the rate Thus the convergence is typically slow. In this 

section, we provide another technique to estimate the distance to points which we 
find, for a given nonzero h, tends to provide better distance estimates. 

The idea of this alternative distance estimation method is to once again check 
for overlap between the sets capp(j/, s, h) and capp( 2 , s, r), but instead of seeking the 
minimum r for which these wave caps overlap, we seek r for which Vol^(capp(y, s, h)n 
capp( 2 ;, s, r)) is half of Vol^(capp(j/, s, h)). 

Before proving that our alternative distance estimation procedure is valid, we 
provide a lemma that shows that the diameter of a wave cap vanishes as the height 
of the cap goes to zero. 

Lemma 4.5. Let y,z gT, s G {0,ar{y))- Then, 

(4.10) lim diam(cap( 2 /, s, hj) = 0. 

h^O 


Proof. Suppose the claim were false. Then there exists a sequence of positive 
real numbers f 0 and points Pi G cap-p{y, s, hi) such that d{x{y, s),Pi) 0. Since 
M is compact, the sequence {pt} has a convergent subsequence. Relabeling this 
subsequence by pt we have that there exists p G M such that pi —> p. But this implies 
that d(p, x{y, s)) ^ 0, hence p ^ x{y, s). On the other hand, since pi G capp(y, s, hi) 
we must have that p G (^^,> 0 ■Sj ^)j but this gives a contradiction, since by 
Lemma 4.2 this implies that p = x(y, s). □ 

We now present our alternative distance estimation method. 

Lemma 4.6. Let y,z G F, s G (0,(Tr(?/)), and 0 < h < crr(y) — s. Let be the 
solution to, 


(4.11) Vol;,(capr( 2 /,s,/i) ncapr(z,s,rft)) = ^ Vol^(capr(?/, s, h)). 


Then, for dh '■= s + r^, we have that dh —t d{z, x(y, s)) as h ^ 0. 

Proof. First, we recall that for s and h as above, capp( 2 /, s, h) will contain a 
non-empty open set, hence the right-hand side of ( |4.11 1 will be nonzero. Thus, from 
the definition of Vh, we conclude that capp(?/, s, h) O capp( 2 ;, s, r^i) is a non-empty 
and proper subset of capp(j/, s, h). Using the definition of capp(z, s, r/j) and r^ we 
conclude that s + ru > d(z, capp(y, s, h)). On the other hand, since the intersection 
between the wave caps is a proper subset of capp(j/, s, h) we see that there exists 
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p C capp(j/, s,/i) \ capp(z, s, r/i). In particular, this implies that s + Vh < d{z,p) < 
dist(z, capp(y, s, h)) + diam(capp(?/, s, h)). Hence, 

d{z, capp(j/, s, h)) < dfi < d{z, capp(j/, s, h)) + diam(capp( 2 /, s, h)). 

Since (i(z, capp(i/, s, ft.)) —)■ d{z,x{y, s)) and diam(capp(j/, s, ft)) —)■ 0 as ft —)■ 0, we 
conclude that dh —t d{z, x{y, s)) as ft —)■ 0. □ 

We summarize the steps of the proof in an algorithmic form in Algorithm [l] 

let: y, z GT and s > 0 with r^t^ys^^z) < T. 

let: fto > 0 small enough that s + ft-o < min{(Tp(?/), T}. 

for 0 < ft < ftg: 

for 0 < r < T — s: 

let: Ti = sip, T 2 = T®+^, Tg = r 4 = Ti V r 2 V Tg 

for a > 0: 

for i = 1,..., 4: 
let: fa,i solve: 

{Kr^+a)PrJ = Pr,b 

for i = 1,... ,4: 
compute: 

TTliri) = 

compute: 


^target cap(ft) := m(T2)-m(Ti) 

Woveriap(ft, r) := m{Ti) - m^Ts) - m(T2) + m{Ti) 

compute: by either: 

method 1: let rt satisfy: 


Th = inf{r > 0 : moveriap(^, ?') > 0 }. 


method 2: let Xh be the solution to: 


compute: r^(y^s){z) by: 


I^overlap (fti I’) — 2 cap (ft) ■ 


rx{v,s){z) = s + 1 ™ ?'/»• 

Algorithm 1: Continuum level description of distance determination algorithm. 


5. Computational experiment. In this section we present a numerical exam¬ 
ple that demonstrates the distance determination procedure that we have described 
in the previous sections. For computational simplicity, we demonstrate our procedure 
in the Euclidean setting. However, we stress that our method can be applied in the 
general Riemannian setting. 
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5.1. Numerical method for the direct problem. For our numerical exam¬ 
ple, we take the manifold M to be the 2-dimensional Euclidean lower half-space 
equipped with the canonical metric and a unit weight function, fi = 1. Under these 
particular choices, the weighted Laplace-Beltrami operator reduces to the Euclidean 
2-dimensional Laplacian. Hence, for our example, the Riemannian wave equation 
( 2 . 1 ) simplifies to the standard 2 + 1 -dimensional wave equation with constant sound- 
speed, c = 1. In order to simulate the situation of partial, local illumination, for our 
source/receiver set, F, we take F = [—L,L] x {0} C DM with L = 2.232. We simulate 
waves propagating for 2T time units, where T = 1.249. 

For sources, we use a basis of Gaussian pulses with the form 


= C exp (-at(< - tif - a^{x - Xjf) , 

with parameters Ot = = 4 • 10^, and where we have selected the constant C to 

normalize the Sources are applied on the regular grid: 


^ 1 0 )) ■ 


is,i — G,i 4“ 1)^G 

Xsj = Xs,i + (j - l)Aa:s 


i = 

J = 1, . . . , 


where the source offset Axs and time between source applications AG are both se¬ 
lected as Axs = Atg = .0147. At each of the = 309 source positions we apply 
Nt,g = 78 sources. For each basis function, we record the Dirichlet trace data on the 
regular grid: 


/. / rv\\ . — G,i 4” l)At,. I — 1,..., 

t r,/) \Xr,k, U • fc = 1, . . . , 

The receiver offset Axg has been taken to be half the source offset, resulting in 
Nx,r = 633 receiver positions. The time between receiver measurements, Atg, is 
1/10 the source time between source applications, resulting in Nt^g = 1701 receiver 
measurements at each receiver position. 

We discretize the Neumann-to-Dirichlet map by solving the forward problem for 
each source and recording its Dirichlet trace at the receiver positions and times 
described above. That is, we compute the following data, 


(5.1) 


Xg^k) — (G,U ■ 


^ — 1 ) • ■ • ; s; j — I 5 ■ • ■ 5 A 3 , s; 1 
1 = 1,..., Nt,g, k= 1,..., Ng,^g j 


To perform the forward modelling, we use the continuous Galerkin finite element 
method with piecewise linear Lagrange polynomial elements and implicit Newmark 
time-stepping. In particular, we use the FEniCS package [IS]. We use a regular 
triangular mesh, where the time step and mesh spacing are selected so that 8 points 
per wavelength (in directions parallel to the grid axes) are used at the frequency fo 
where the spectrum of the temporal portion of the source falls below 10 “® times its 
maximum value. 


5.2. Solving the control problem. We discretize the connecting operator K 
by approximating its action as an operator on span{(^ij }. That is, we use the discrete 
Neumann-to-Dirichlet data, (5.1), to discretize Kg by formula (|3.2|), where t = T. 


To be specific, we first compute the Gram matrix [G]ij = ((pi,ipj) and its inverse 
and then compute the matrices for the operators JA^, RAf and RJ acting 
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(a) 


(b) 


(c) 


(d) 


(e) 



A 



^h400 

200 


20 

10 

0 


0.5 

0 


1.5 


0.5 


0.6 

0.4 



Table 5.1: Demonstration of the action of the operators appearing in (5.2), (a) a 
source /, (b) RJf, (c) RAfRJf, (d) JAp^/, and (e) Kf. Note that all plots are 
on the function level, and each function depicted was obtained from coefficients that 
were computed as described in Section 15.21 


on span{(/?i j} as follows: 


[JAfjij = %kiVk,JAf'(pj), 

k 

[i?Ap]ij = y^[G ^]ik{‘Pk, RA^ipj), 

k 


Using these operators we can compute the matrix for AT, 
(5.2) [K] = [JAl^]-[RAl][RJ]. 


In Table 5.1 we demonstrate the effect of the constituent matrices in (5.2) on a source 


/ and how these combine to form Kf. 

For r € G(r), with 0 < r < T, we obtain the matrix [Kt] discretizing the con¬ 
necting operator Kt by masking the entries in [K\ that correspond to basis functions 
ipi^j with centers (ts.u Xsj) ^ Sr- We note that, in practice, we find that this tends to 
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provide a better approximation to Kt than computing the matrix [Pr] and computing 
the product 

We consider the discretized control problem 


(5.3) 


i[Kr]+a[Pr])[U] = [Pr][b], 


where we use the matrix [Pt] to refer to the mask described above, and use a = 10“®. 
Recall that b is the function given by b{t,x) — T — t, as defined beneath (3.81. To 
solve (5.3) for [fa], we use restarted GMRES. In Figure and Table 5.2 we depict 
control solutions fa = their associated wavefields u^’=‘{T, •). A volume 

estimate 7ti(t) for m(T) is obtained from [fa] by comput ing th e discretized inner 
product TO(r) = which approximates m{T) as in (3.10). For the remainder 

of this paper we will continue to use the notation TO(r) to indicate the approximation 
to mir) computed in this fashion. 



Fig. 5: Illustration of a source fa and the corresponding wa vefie ld (T, •) for which 

Ilf ‘ 


‘(T, •) « 1m(t)- Here, r corresponds to T 4 ^i from Table 


5.2 


both plots with the same horizontal axis, we have extended" 
[o,r] xT. 


Note that, to show 
to zero outside of 


5.3. Estimating distances. We estimate distances between z G T and points 
of the form x{y,s) where y = (0,0). In particular, for each fixed s we estimate the 
distances d{x{y, s), {zi, 0)) for uniformly spaced [zi, 0) G [—1.175,1.175] x {0} C T. We 
take the offset Az between the points Zi equal to Az = dAa;^ = .0588, and select the 
points {zi,0) to coincide with every fourth source position. As a proxy for estimating 
the distance to x{y, s), we use a target wave cap of the form capp(j/, s, h) with height 
h = .025, and estimate the distances for s = .125, .25, .375, .5. 


For each s we solve the discrete control problem (5.3) in order to obtain estimates 
TO(slr) and fh{Ty~^^) for the respective volumes of M(r,s) and M{y,s + h). From 
these, we estimate the volume of the target cap by, 

TOtarget cap = Ai(t®+'‘) - m(slr). 

For each point (zi,0), we also solve control problems to obtain volume estimates 
TO(r^’^ 0 )) where we select the parameters rj,j = 1,... ,Nr 

so that the two sets {s + rj : j = 1,..., W} = : ts,k > r} coincide. We implement 

the distance estimation procedure described in Lemma 4.6 to estimate r„,(y^s-)((zi,0)) 






















CONSTRUCTING TRAVEL TIME DISTANCES 


17 



Table 5.2: Illustration of the essential features of sources and waves used in distance 
estimation procedure. Solutions, fa, to the discretized control problem plotted next 
to their associated wavefields approximating Plotted for r of the 

form Ti = sir, T 2 = and T 4 j = ti V T 2 V ts^, where y = (0,0), 

z = (0.529, 0),h = .05, and = 0.118,0.235,0.338. 


as follows: for each rj we estimate the volume of capr(?/, s, h) n capr((zi, 0), s,rj) by 
first computing. 


^overlapj = V V slp) - - m(Ty +*) + 77l(slr). 
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Table 5.3: Plots of relative overlap volumes, uioveriap /bi. target cap, vs. r. Plotted for 
s = .25, h = .025 and (clockwise from the top left) z = 0.176 ■ i, i = 0,... ,3. The 
markers denote the relative overlap volumes estimated in the distance determination 
procedure, and the lines indicate the analytical relative overlap volumes. 


We then find the indices j, j + 1 for which 


(5.4) 


hioverlap,j ^ 


1 , 

2 ^target cap ^ ^overlap,j-t- 1 , 


and estimate by linearly interpolating between Vj and rj^i. This procedure ap¬ 
proximates (4.11). We depict the results of the volume overlap estimation in Table 
|5.3[ Since the volumes in these images have all been normalized by the target cap 
volumes, computing rh by (|5.4|) corresponds to hnding the x-value where the curve 


connecting the data points passes through the line y = 0.5. We depict the distance 
estimation results in Figure]^ 


5.4. Discussion of sources of numerical errors and instability. Examining 
Figure one can see that in each of the estimated distance curves, the distances 
are over-estimated for z = (zi,0) near y = (0,0). This error results in part from 
the distance estimation method. For example, when z = y the correct distance 
d = r + s would be obtained by taking r = 0. On the other hand, when z = y, 
both of the wave caps used in the distance estimation procedure are centered on the 
same point, so for 0 < r < h the variable wave cap, capp(z, s,r), coincides with 
capp(z, s,r) n capp( 2 /, s,/i). From the definition of r^, we find that we will have 
0 < rh < h. Thus the distance estimate dh will necessarily over-estimate d{y, x{y, s)). 
Similar remarks apply for estimating d((zi, 0), a;(j/, s)) for (zi,0) near y, although the 
strength of this effect decreases as (zi,0) gets further from y. We call this source of 
error geometric distortion, since it results entirely from the geometry of our distance 
estimation procedure and is independent of errors arising from the control problems. 
In Figure we depict the geometric distortion by repeating our distance estimation 
technique with exact volume measurements. Note that the distances in Figure are 
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Fig. 6: Distance estimates (markers) for d{x{y,s),{zi,0)) for y = (0,0) and s = 
.125, .25, .375, .5, plotted along with the true distances (solid curves). 


overestimated at all points, which contrasts most with the distances estimated at large 
offsets in Figure]^ 



Fig. 7: Demonstration of geometric distortion. Distances (markers) are esti¬ 
mated by using the distance estimation technique on exact volumes and plotted for 
d{x{y, s), {zi, 0)) for y = (0,0) and s = .125, .25, .375, .5, along with the true distances 
(solid curves). 
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Our numerical tests suggest that the dominant source of error comes from the 
control step. In order to discuss this instability, we return to considering the con¬ 
tinuum problem. Taking r G C(r), we can ask whether there exists / S H^{St) for 
some s S K for which Wrf = This question can be answered by considering 

the more general problem of exact controllability, in which one seeks to determine 
when the equation = {wq,wi) has a solution in H^{St) for any 

{wq,Wi) belonging to an appropriate space of Cauchy data for the wave equation. 

In jl], the question of exact controllability is considered. One of the main results 
of that paper is that the ray geometry of the wave equation can be used to determine 
necessary and sufficient conditions for exact controllability. Using the same set of 
ideas to those in [1], it is shown in [5] that in order for exact controllability to hold 
for Wr in M[t) from Sr, the following geometric controllability condition must hold: 

Each generalized bicharacteristic {x{t),t) satisfying x(T) £ M(t), passes over 
Sr id S'r in a non-diffractive point. 

Here, S'r = {{t,x) G T x {T,2T) : T < t < T t{x)}. We recall that Sr is 
defined by (2.3), and note that S'r is the temporal reflection of Sr across t = T. For 
a generalized bicharacteristic {x{t),t), the path x{f) is a unit speed geodesic in the 
interior of M and it is reflected according to Snell’s law when it intersects the boundary 
dM transversally. Tangential intersections with the boundary can cause the path to 
glide along the boundary, and in the case of an infinite-order contact, the path x{t) 
can be continued in many ways, see [1]. We refer also to [2] for the definition of 
non-diffractive points. The geometric controllability condition is necessary for exact 
control to hold from Sr, since when it fails for {x,£f) G S*{M{t)), propagation of 
singularities implies that for any s G K and any / G H'^{Sr), {x,f,) ^ WF(u'^(T, •)), 
see e.g. [TSl Section 23]. Here, WF(u'^(r, •)) denotes the wave front set of 
and we refer to na Def. 8.1.2] for its definition. We have also provided the definition 
of wave front set in Appendix 0 Thus, if rc G L^(M(t)) has {x,f) G WF(ri;) then, 
for each s G K, there does not exist / G H'^{Sr) for which Wrf = w. 

In our computational experiment, the geometric controllability condition actually 
fails over every point in Mir). This is due to the fact that at each x G M{t) there 
exists a family of unit-speed geodesic rays with = {x,f,) G Sr{M) and 

^ belonging to a cone over x, for which the corresponding geodesics 7 fail to pass 
over S'.r U S'].. In our computational experiment, we observe instabilities near those 
X G M(t) where WF(1 m(t)) meets the cone over which exact control fails. In the 
case of r = V sir, these effects occur where 9 M(t) fails to be (7°° smooth and 
{a:} X (M" \ 0) C WF(lM(r)) ■ We refer to Appendix [A| for further analysis. In 
particular these effects occur for x in the bottom left and right edges of capr(y, s, h), 
where dMir) fails even to be . In addition, we similarly observe instabilities near 
the points (±T, —s), where the flat portion of dMir) transitions into a circle and fails 
to be C2. _ 

by plotting a wavefield u^(T,-) ap- 


We demonstrate these effects in Figure 


proximating 1m(t) for y = (0,0), s = .5, and h = .2. The former instabilities occur 
near the points (±.5,—.5), and the latter instabilities occur near (±1.15,—.5). We 
contrast this with the case where r = which we show in Figure 


8 b 


In this 


second example, the domain of influence is a disk and every co-vector in Wi''(lM(T)) 
can be controlled, and unlike the first example, we observe no instabilities. Note that 
in all of the examples in Figure we use a smaller F than in our distance calculations, 
using L = 1.153 and a finer basis with ot = Oa, = 16 • 10^. 

we plot the wavefield u^iT, •) that approximates iM(sir)- ^^ote that 


In Figure 


8 c 
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Fig. 8: (a) wavefield demonstrating instability of the solution to the control problem 
when WF(lM(r)) contains uncontrollable directions over (±.5,—.5) and (±1.15,—.5) 
(b) A wavefield for which all directions in WF(ljv^(,-)) are controlled, (c) Another 
waveheld demonstrating instability, with uncontrollable directions in WF(ljv^(^)) over 
(±1.15,—.5). (d) The difference between the wavefields in (a) and (c), note that this 
corresponds to an approximation to ^M(ca.-p(y,s,h)) as used in the distance estimation 
procedure. Moreover, the instabilities in (a) and (c) located over (±1.15,—.5) cancel 
each other. 


as in the case of r = V sir, we observe instabilities near the points (±1.15, —.5). 
In Figure [8dj we plot the difference between the wave helds approximating lM(T“+'‘vsir) 
and lM(sir)j and note that this difference yields an approximation to the characteristic 
function of capp(j/, s, h). In particular, notice that the instabilities observed near 
(±1.15, —.5) in Figures 8a and 8c completely cancel in Figure 8d Since our distance 


determination relies primarily on the volumes of wave caps, which are obtained by 
taking differences in this fashion, we find that the instabilities near the cap bases tend 
to provide the main source of error for our distance estimation procedure. 


6. Conclusions. In this paper we have demonstrated a method to construct 
distances between boundary points and interior points with fixed semi-geodesic co¬ 
ordinates. The procedure is local in that it utilizes the local Neumann-to-Dirichlet 
map for an acoustic wave equation on a Riemannian manifold with boundary. Our 
procedure differs from earlier results in that it utilizes volume computations derived 
from local data in order to construct distances. Finally, we have provided a computa- 
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tional experiment to demonstrate the implementation of our distance determination 
procedure and we have discussed the main sources of numerical error. 

Appendix. Wave front set of 

In Section all of the functions r that we consider give rise to sets M (r) with 
piecewise smooth boundary. For all such t, if x G dMir) is a point for which dMlr) 
is not smooth at x, then either dM{T) is not at x, or dM{T) fails to be at x. 
In this section, we examine the former case, by computing the wave front set of Im(t) 
over a point x where dM{T) fails to be . The other case is similar and we omit it. 
In particular, we will show the following lemma: 

Lemma A.l. Let h : M. ^ M. be continuous, and suppose that h is smooth on 
(—00,0) and on (0, oo), and that h does not belong to Let A be the set, 

A = {(xi,x2) G K 2 . ^2 < Then, {0} x (K^ \ q) C WF(U). 

From the lemma, it follows that if r is one of the functions considered in Section 
and dM(r) fails to be at x G 9M(r), then {x} x \ 0) C WF(1m(t))- Put 
differently, WF(Ijv/(t)) contains all cotangent directions above the point x. 

Before proceeding with the proof, we recall the definition of wave front set, 

Definition A. 2. Let X c i?" open. If u G DfX), then the closed subset of 
X X (M" \ 0) defined by, 

WF(u) = {(x, 0 G A X (K" \ 0) : C e S,(u)} 
is called the wave front set of u, where £ 2 ,( 11 ) is the set 

£,(u) = nS(0u), </> G C'o“(A), cf{x) ^ 0, 

0 


and, for v G 'DfX), £(u) is the complement of the set of p in K" \ 0 for which there 
exists a conic neighborhood V of rj such that for each A G N there exist Cjv > 0 such 
that, for all ^ G V, 




To prove Lemma |A.l[ we require the following result. 

Lemma A.3. Let h G M., u G C'(j“(K). Suppose that </> G C°°(M) is real valued and 
has no critical points in supp(i(). Then 



dx = 


\4>'{h) 


X'^cffh) A3 


where |i?| is bounded by a constant that depends only on supp(it), min 2 ;gsupp(ti) l<(’^(a;)| 
and the norms of u and 4>. 

Proof. We define the differential operator L = i(j)'~^dx. Then 



M' 



uLe *^‘^dx 


+ 


I [ {L*u)e-^^^dx. 

^ J — 00 


x—h 
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Notice, L*u = —idx{u/(j)'). Then, 


1 

A 


f {L*u)e-^^^dx = ^ f {L*u)Le-’-^‘<’dz 

J —oo ^ J —oo 


dx{u/(j)')e 






c—h 


A 2 


{{L*fu)e-^^^dx. 


Finally, 


^ £ {{L*fu)e-^^Ux = ^ £ {{L*fu)Le-^^^dz 




' —OO 
ph 


A3 

x—h J-oo 


{{L*fu)e-^^'^dx. 


Proof, [of Lemma A.l Let u G and suppose that u = 1 near the origin. 

We consider the Fourier transform 


mIa(A^) = f u{x)lA{x)e ’‘^^^dx 
dR 2 


/ +00 ^h(x^) 

/ u{x)e-^^^^^"dx^dx\ 

-OO ^ —OO 


where ^ is a unit vector and A > 0. 

Suppose first that ^2 0. Then 

/ u{x)e dx^ = •—^—— -- 

J —00 ^^2 


A 2 ^i 


where w{x^) = u{x^,h{x^)) and v{x^) = dx 2 u{x^, h{x^)). Note that R is compactly 
supported with respect to x^ since u is. Note also that w and v are compactly 
supported and w = 1 near the origin. If supp(m) is small then supp(w) and supp(u) 
are also small. Moreover, if supp(r(;) U supp(?;) is small enough then h is smooth in 
(supp(w) U supp(u)) \ 0. In particular, w and v are smooth then. 

We define ^(x^) = fix^ +^ 2 h(x^), r/>± = 4‘\±x^>o^ and define also h± analogously. 
Suppose that (/>'_(0) = 0. Then </> has no critical points in the set 

{x^ < 0} n supp(w) if supp(rc) is small enough. Therefore 


we-^^^dx^ = 


Likewise, if also </>Y(0) ^ Oj then 


^+00 


A</>'_(0) 


— I 


we-^^^^dx^ = 

/o A0'+(O) 


i?+(A, h, (f), u) 
A 2 ■ 


R- (A, h, cj), u) 
A 2 • 


So, 


^ + 00 


we = i 




X-^ + 0{X-^). 
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An analogous argument applies to u, showing that, 



O(A-i). 


Hence, we conclude, 

Thus does not decay rapidly if </)+(0) ^ ^^'^-(0) which again is equivalent to 

To summarize, if /i+(0) 7 ^ h'_{0) then all the directions except possibly (1,0) and 
the four directions 

(-/l±(0)C^?^), where |^ 2 | = (|^±(0)P + 1)"^/^, 

are in I]o(1a). Finally, since Eo(lyi), is a closed conic subset of \ 0, we conclude 
that Eo(1a) = \ 0, and hence {0} x (M^ \ 0) C WF(l^). □ 
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